%% bdary_node (for Dirichlet)
function bdn = bdary_node(mesh, fespace, index)
bde = bdary_edge(mesh, index);
bdn = unique(bde(:,2:3), "stable");
if 2 == size(bdn,2)
    bdn = reshape(bdn, [], 1);
end
if ismember(fespace.basis_type, ["P2", "P3", "P4", "P5"])
    basis_type = char(fespace.basis_type);
    r = str2double(basis_type(2))-1;
    node_ed = zeros(r, size(bde,1));
    for k = 1:r
        node_ed(k,:) = bde(:,4)'+mesh.Np+(k-1)*mesh.Ned;
    end
    if bde(1,2) == bde(end,3)
        bdn = reshape([bdn'; node_ed], [], 1);
    else
        O = zeros(r, 1);
        bdn = reshape([bdn'; [node_ed, O]], [], 1);
        bdn = bdn(1:end-r,:);
    end
end
end